"""
# Parameter Optimization and Research Workflows

This module provides comprehensive tools for parameter sweeps, optimization, 
sensitivity analysis, and reproducible research workflows for PDEJulia.jl.
It's designed to accelerate research by automating parameter exploration
and enabling systematic studies.

## Key Features

- Multi-dimensional parameter sweeps with parallel execution
- Bayesian optimization for expensive simulations  
- Sensitivity analysis and uncertainty quantification
- Automatic result storage and retrieval
- Publication-ready visualization and reporting
- Integration with common PDE research workflows

## Optimization Methods

- Grid search for exhaustive exploration
- Random sampling for high-dimensional spaces
- Latin hypercube sampling for space-filling designs
- Bayesian optimization with Gaussian processes
- Multi-objective optimization (Pareto fronts)
- Gradient-based optimization where applicable

## Applications

- Parameter fitting to experimental data
- Optimal control of PDE systems
- Design optimization (materials, devices)
- Sensitivity and uncertainty analysis
- Model selection and comparison
- Reproducible research workflows

## References

- Jones, D. R. (1998). Efficient global optimization of expensive black-box functions.
- Sobol, I. M. (2001). Global sensitivity indices for nonlinear mathematical models.
- McKay, M. D. (1979). Comparison of three methods for selecting values of input variables.
"""

using LinearAlgebra, Statistics, Random, Printf
using Distributed, SharedArrays  # For parallel computation
using JSON, HDF5, CSV, DataFrames  # For data management
using Plots, StatsPlots  # For visualization
using Optim  # For optimization algorithms

export ParameterSpace, OptimizationConfig, ResearchWorkflow
export parameter_sweep, bayesian_optimization, sensitivity_analysis
export create_research_workflow, reproduce_study, generate_report

"""
    ParameterSpace

Defines a multi-dimensional parameter space for exploration.

# Fields
- `parameters::Dict{Symbol, Any}`: Parameter definitions
- `constraints::Vector{Function}`: Parameter constraints  
- `fixed_params::Dict{Symbol, Any}`: Fixed parameters
- `metadata::Dict{String, Any}`: Additional information
"""
struct ParameterSpace
    parameters::Dict{Symbol, Any}
    constraints::Vector{Function}
    fixed_params::Dict{Symbol, Any}
    metadata::Dict{String, Any}
end

"""
    OptimizationConfig

Configuration for optimization algorithms.

# Fields
- `method::Symbol`: Optimization method
- `objective::Function`: Objective function to minimize
- `max_iterations::Int`: Maximum iterations/evaluations
- `tolerance::Float64`: Convergence tolerance
- `parallel::Bool`: Use parallel evaluation
- `options::Dict`: Method-specific options
"""
struct OptimizationConfig
    method::Symbol
    objective::Function
    max_iterations::Int
    tolerance::Float64
    parallel::Bool
    options::Dict{String, Any}
end

"""
    ResearchWorkflow

Encapsulates a complete research workflow with reproducibility.

# Fields
- `name::String`: Workflow identifier
- `description::String`: Detailed description
- `parameter_space::ParameterSpace`: Parameter definitions
- `simulation_function::Function`: Main simulation function
- `analysis_functions::Vector{Function}`: Post-processing functions
- `storage_config::Dict`: Data storage configuration
"""
struct ResearchWorkflow
    name::String
    description::String
    parameter_space::ParameterSpace
    simulation_function::Function
    analysis_functions::Vector{Function}
    storage_config::Dict{String, Any}
end

"""
    ParameterSpace(parameters::Dict; constraints=Function[], fixed_params=Dict(), metadata=Dict())

Create parameter space definition.

# Examples
```julia
# Define parameter ranges
params = Dict(
    :diffusion => 0.01:0.01:0.1,
    :reaction_rate => [0.1, 0.5, 1.0, 2.0],
    :grid_size => [64, 128, 256],
    :epsilon => LogRange(0.001, 0.1, 10)
)

# Add constraints
constraints = [
    p -> p.diffusion * p.reaction_rate < 1.0,  # Stability constraint
    p -> p.grid_size >= 64                      # Resolution constraint
]

space = ParameterSpace(params, constraints=constraints)
```
"""
function ParameterSpace(parameters::Dict; constraints=Function[], 
                       fixed_params=Dict(), metadata=Dict())
    return ParameterSpace(parameters, constraints, fixed_params, metadata)
end

"""
    OptimizationConfig(; method=:bayesian, objective=nothing, max_iterations=100,
                       tolerance=1e-6, parallel=true, kwargs...)

Create optimization configuration.
"""
function OptimizationConfig(; method=:bayesian, objective=nothing, max_iterations=100,
                           tolerance=1e-6, parallel=true, kwargs...)
    
    if objective === nothing
        throw(ArgumentError("Objective function must be specified"))
    end
    
    options = Dict(string(k) => v for (k, v) in kwargs)
    return OptimizationConfig(method, objective, max_iterations, tolerance, parallel, options)
end

"""
    generate_parameter_combinations(space::ParameterSpace; method=:grid, n_samples=nothing)

Generate parameter combinations for exploration.

# Methods
- `:grid`: Full factorial grid (exhaustive)
- `:random`: Random sampling
- `:lhs`: Latin hypercube sampling
- `:sobol`: Sobol sequence (quasi-random)
"""
function generate_parameter_combinations(space::ParameterSpace; method=:grid, n_samples=nothing)
    if method == :grid
        return generate_grid_combinations(space)
    elseif method == :random
        n_samples = n_samples === nothing ? 1000 : n_samples
        return generate_random_combinations(space, n_samples)
    elseif method == :lhs
        n_samples = n_samples === nothing ? 100 : n_samples
        return generate_lhs_combinations(space, n_samples)
    elseif method == :sobol
        n_samples = n_samples === nothing ? 1024 : n_samples
        return generate_sobol_combinations(space, n_samples)
    else
        throw(ArgumentError("Unknown sampling method: $method"))
    end
end

"""
    generate_grid_combinations(space::ParameterSpace)

Generate full factorial parameter grid.
"""
function generate_grid_combinations(space::ParameterSpace)
    param_names = collect(keys(space.parameters))
    param_values = [collect(values) for values in values(space.parameters)]
    
    combinations = []
    
    for combo in Iterators.product(param_values...)
        param_dict = Dict(zip(param_names, combo))
        
        # Merge with fixed parameters
        full_params = merge(space.fixed_params, param_dict)
        
        # Check constraints
        if all(constraint(full_params) for constraint in space.constraints)
            push!(combinations, full_params)
        end
    end
    
    return combinations
end

"""
    generate_random_combinations(space::ParameterSpace, n_samples::Int)

Generate random parameter combinations.
"""
function generate_random_combinations(space::ParameterSpace, n_samples::Int)
    combinations = []
    param_names = collect(keys(space.parameters))
    
    for _ in 1:n_samples * 2  # Generate extra to account for constraint rejections
        param_dict = Dict{Symbol, Any}()
        
        for (name, values) in space.parameters
            if isa(values, AbstractRange) || isa(values, AbstractVector)
                param_dict[name] = rand(values)
            else
                # Assume it's a distribution or callable
                param_dict[name] = rand(values)
            end
        end
        
        # Merge with fixed parameters
        full_params = merge(space.fixed_params, param_dict)
        
        # Check constraints
        if all(constraint(full_params) for constraint in space.constraints)
            push!(combinations, full_params)
            if length(combinations) >= n_samples
                break
            end
        end
    end
    
    return combinations[1:min(n_samples, length(combinations))]
end

"""
    generate_lhs_combinations(space::ParameterSpace, n_samples::Int)

Generate Latin hypercube sampling combinations.
"""
function generate_lhs_combinations(space::ParameterSpace, n_samples::Int)
    param_names = collect(keys(space.parameters))
    n_params = length(param_names)
    
    # Generate LHS samples in [0,1]^n
    lhs_samples = latin_hypercube_sample(n_samples, n_params)
    
    combinations = []
    
    for i in 1:n_samples
        param_dict = Dict{Symbol, Any}()
        
        for (j, (name, values)) in enumerate(space.parameters)
            # Map [0,1] to parameter range
            u = lhs_samples[i, j]
            
            if isa(values, AbstractRange)
                # Linear interpolation for ranges
                idx = 1 + floor(Int, u * (length(values) - 1))
                idx = clamp(idx, 1, length(values))
                param_dict[name] = values[idx]
            elseif isa(values, AbstractVector)
                # Discrete selection
                idx = 1 + floor(Int, u * length(values))
                idx = clamp(idx, 1, length(values))
                param_dict[name] = values[idx]
            else
                # Assume continuous distribution
                param_dict[name] = quantile(values, u)
            end
        end
        
        # Merge with fixed parameters
        full_params = merge(space.fixed_params, param_dict)
        
        # Check constraints
        if all(constraint(full_params) for constraint in space.constraints)
            push!(combinations, full_params)
        end
    end
    
    return combinations
end

"""
    latin_hypercube_sample(n_samples::Int, n_dims::Int)

Generate Latin hypercube samples in [0,1]^n.
"""
function latin_hypercube_sample(n_samples::Int, n_dims::Int)
    samples = zeros(n_samples, n_dims)
    
    for j in 1:n_dims
        # Generate random permutation
        perm = randperm(n_samples)
        
        # Generate stratified samples
        for i in 1:n_samples
            samples[i, j] = (perm[i] - rand()) / n_samples
        end
    end
    
    return samples
end

"""
    parameter_sweep(space::ParameterSpace, simulation_function::Function;
                   method=:grid, n_samples=nothing, parallel=true, 
                   storage_file=nothing, verbose=true)

Perform comprehensive parameter sweep.

# Arguments
- `space`: Parameter space definition
- `simulation_function`: Function to evaluate for each parameter set
- `method`: Sampling method (:grid, :random, :lhs, :sobol)
- `n_samples`: Number of samples (for non-grid methods)
- `parallel`: Use parallel computation
- `storage_file`: File to store results (HDF5 or JSON)
- `verbose`: Print progress information

# Returns
- `NamedTuple` with parameter combinations and results
"""
function parameter_sweep(space::ParameterSpace, simulation_function::Function;
                        method=:grid, n_samples=nothing, parallel=true,
                        storage_file=nothing, verbose=true)
    
    if verbose
        println("🔍 Parameter Space Exploration")
        println("=============================")
        println("Method: $method")
        println("Parallel: $parallel")
    end
    
    # Generate parameter combinations
    combinations = generate_parameter_combinations(space, method=method, n_samples=n_samples)
    n_combinations = length(combinations)
    
    if verbose
        println("Generated $(n_combinations) parameter combinations")
        println("Estimated time per simulation: depends on problem size")
        println()
    end
    
    # Storage for results
    results = Vector{Any}(undef, n_combinations)
    errors = Vector{Union{Nothing, Exception}}(undef, n_combinations)
    
    if parallel && nprocs() > 1
        # Parallel execution
        if verbose
            println("Starting parallel execution on $(nprocs()) processes...")
        end
        
        @sync @distributed for i in 1:n_combinations
            try
                result = simulation_function(combinations[i])
                results[i] = result
                errors[i] = nothing
            catch e
                results[i] = nothing
                errors[i] = e
                if verbose
                    println("Error in combination $i: $e")
                end
            end
        end
    else
        # Serial execution with progress
        if verbose
            print("Progress: ")
        end
        
        for (i, params) in enumerate(combinations)
            try
                results[i] = simulation_function(params)
                errors[i] = nothing
            catch e
                results[i] = nothing
                errors[i] = e
                if verbose
                    println("\nError in combination $i: $e")
                end
            end
            
            if verbose && i % max(1, n_combinations ÷ 20) == 0
                progress = i / n_combinations * 100
                print("$(Int(round(progress)))% ")
            end
        end
        
        if verbose
            println()
        end
    end
    
    # Filter successful results
    successful_mask = errors .=== nothing
    successful_combinations = combinations[successful_mask]
    successful_results = results[successful_mask]
    n_successful = sum(successful_mask)
    
    if verbose
        println("✅ Parameter sweep completed!")
        println("Successful runs: $n_successful / $n_combinations")
        if n_successful < n_combinations
            n_failed = n_combinations - n_successful
            println("Failed runs: $n_failed ($(round(n_failed/n_combinations*100, digits=1))%)")
        end
        println()
    end
    
    # Create comprehensive results structure
    sweep_results = (
        parameters = successful_combinations,
        results = successful_results,
        errors = errors,
        metadata = Dict(
            "method" => method,
            "n_total" => n_combinations,
            "n_successful" => n_successful,
            "parameter_space" => space,
            "timestamp" => string(now())
        )
    )
    
    # Store results if requested
    if storage_file !== nothing
        save_sweep_results(sweep_results, storage_file, verbose)
    end
    
    return sweep_results
end

"""
    save_sweep_results(results, filename::String, verbose=true)

Save parameter sweep results to file.
"""
function save_sweep_results(results, filename::String, verbose=true)
    ext = lowercase(splitext(filename)[2])
    
    if ext == ".h5" || ext == ".hdf5"
        save_results_hdf5(results, filename)
    elseif ext == ".json"
        save_results_json(results, filename)
    else
        throw(ArgumentError("Unsupported file format: $ext"))
    end
    
    if verbose
        println("Results saved to: $filename")
    end
end

"""
    save_results_hdf5(results, filename::String)

Save results in HDF5 format.
"""
function save_results_hdf5(results, filename::String)
    h5open(filename, "w") do file
        # Save metadata
        metadata_group = create_group(file, "metadata")
        for (key, value) in results.metadata
            if isa(value, String) || isa(value, Number)
                metadata_group[key] = value
            else
                metadata_group[key] = string(value)
            end
        end
        
        # Save parameters (flattened)
        param_group = create_group(file, "parameters")
        for (i, params) in enumerate(results.parameters)
            param_subgroup = create_group(param_group, "run_$i")
            for (key, value) in params
                param_subgroup[string(key)] = value
            end
        end
        
        # Save results (this is problem-specific and may need customization)
        result_group = create_group(file, "results")
        for (i, result) in enumerate(results.results)
            if result !== nothing
                # This is a simplified save - real implementation would handle
                # different result types appropriately
                result_subgroup = create_group(result_group, "run_$i")
                if isa(result, Dict)
                    for (key, value) in result
                        if isa(value, AbstractArray) && eltype(value) <: Number
                            result_subgroup[string(key)] = value
                        elseif isa(value, Number)
                            result_subgroup[string(key)] = value
                        end
                    end
                end
            end
        end
    end
end

"""
    save_results_json(results, filename::String)

Save results in JSON format (simplified).
"""
function save_results_json(results, filename::String)
    # Convert results to JSON-serializable format
    json_data = Dict(
        "metadata" => results.metadata,
        "n_runs" => length(results.parameters),
        "parameters" => [Dict(string(k) => v for (k, v) in params) 
                        for params in results.parameters]
    )
    
    # Note: Full result data may be too complex for JSON
    # This saves only the parameter combinations
    
    open(filename, "w") do file
        JSON.print(file, json_data, 2)
    end
end

"""
    analyze_parameter_effects(sweep_results, response_variable::Symbol; 
                             plot_results=true, n_top=10)

Analyze parameter effects on response variable.

# Arguments
- `sweep_results`: Results from parameter_sweep
- `response_variable`: Symbol identifying the response to analyze
- `plot_results`: Create visualization plots
- `n_top`: Number of top parameter sets to highlight

# Returns
- Statistical analysis of parameter effects
"""
function analyze_parameter_effects(sweep_results, response_variable::Symbol;
                                  plot_results=true, n_top=10)
    
    println("📊 Parameter Effect Analysis")
    println("===========================")
    
    parameters = sweep_results.parameters
    results = sweep_results.results
    
    # Extract response variable values
    response_values = []
    valid_params = []
    
    for (i, result) in enumerate(results)
        if result !== nothing && haskey(result, response_variable)
            push!(response_values, result[response_variable])
            push!(valid_params, parameters[i])
        end
    end
    
    if isempty(response_values)
        println("No valid data found for response variable: $response_variable")
        return nothing
    end
    
    response_values = collect(response_values)
    n_valid = length(response_values)
    
    println("Response variable: $response_variable")
    println("Valid data points: $n_valid")
    println("Response statistics:")
    println("  Mean: $(round(mean(response_values), digits=6))")
    println("  Std:  $(round(std(response_values), digits=6))")
    println("  Range: [$(round(minimum(response_values), digits=4)), $(round(maximum(response_values), digits=4))]")
    println()
    
    # Identify parameter effects
    param_names = collect(keys(valid_params[1]))
    param_effects = Dict{Symbol, Float64}()
    
    for param_name in param_names
        param_values = [p[param_name] for p in valid_params]
        
        # Calculate correlation (for numeric parameters)
        if all(isa(v, Number) for v in param_values)
            correlation = cor(param_values, response_values)
            param_effects[param_name] = abs(correlation)
        end
    end
    
    # Sort parameters by effect magnitude
    sorted_effects = sort(collect(param_effects), by=x->x[2], rev=true)
    
    println("Parameter effects (absolute correlation):")
    for (param, effect) in sorted_effects[1:min(5, length(sorted_effects))]
        println("  $param: $(round(effect, digits=4))")
    end
    println()
    
    # Identify top performing parameter sets
    sorted_indices = sortperm(response_values, rev=true)
    top_indices = sorted_indices[1:min(n_top, length(sorted_indices))]
    
    println("Top $n_top parameter combinations:")
    for (rank, idx) in enumerate(top_indices)
        println("  Rank $rank ($(response_variable) = $(round(response_values[idx], digits=6))):")
        for (param, value) in valid_params[idx]
            println("    $param: $value")
        end
        println()
    end
    
    if plot_results && length(param_names) > 0
        # Create parameter effect plots
        effect_plot = create_parameter_effect_plots(valid_params, response_values, 
                                                   sorted_effects[1:min(4, end)])
        return (effects=param_effects, top_parameters=valid_params[top_indices], 
                response_values=response_values[top_indices], plot=effect_plot)
    end
    
    return (effects=param_effects, top_parameters=valid_params[top_indices],
            response_values=response_values[top_indices])
end

"""
    create_parameter_effect_plots(parameters, responses, top_effects)

Create visualization of parameter effects.
"""
function create_parameter_effect_plots(parameters, responses, top_effects)
    n_effects = length(top_effects)
    plots_list = []
    
    for (param_name, effect) in top_effects
        param_values = [p[param_name] for p in parameters]
        
        if all(isa(v, Number) for v in param_values)
            # Scatter plot for numeric parameters
            p = scatter(param_values, responses, 
                       xlabel=string(param_name), ylabel="Response",
                       title="Effect: $(round(effect, digits=3))",
                       alpha=0.6, ms=3)
            
            # Add trend line
            if length(unique(param_values)) > 1
                plot!(p, param_values, predict_trend(param_values, responses),
                      color=:red, linewidth=2, label="Trend")
            end
            
            push!(plots_list, p)
        else
            # Box plot for categorical parameters
            unique_vals = unique(param_values)
            grouped_responses = [responses[param_values .== val] for val in unique_vals]
            
            p = boxplot(string.(unique_vals), grouped_responses,
                       xlabel=string(param_name), ylabel="Response",
                       title="Effect: $(round(effect, digits=3))")
            
            push!(plots_list, p)
        end
    end
    
    if length(plots_list) > 0
        layout_dims = isqrt(length(plots_list))
        layout_dims = (layout_dims, ceil(Int, length(plots_list)/layout_dims))
        combined_plot = plot(plots_list..., layout=layout_dims, size=(800, 600))
        return combined_plot
    end
    
    return nothing
end

"""
    predict_trend(x_values, y_values)

Simple linear trend prediction.
"""
function predict_trend(x_values, y_values)
    # Simple linear regression
    n = length(x_values)
    x_mean = mean(x_values)
    y_mean = mean(y_values)
    
    slope = sum((x_values .- x_mean) .* (y_values .- y_mean)) / sum((x_values .- x_mean).^2)
    intercept = y_mean - slope * x_mean
    
    return slope * x_values .+ intercept
end

"""
    create_research_workflow(name::String, description::String,
                            parameter_space::ParameterSpace,
                            simulation_function::Function;
                            analysis_functions=Function[],
                            storage_path="./results")

Create a reproducible research workflow.
"""
function create_research_workflow(name::String, description::String,
                                 parameter_space::ParameterSpace,
                                 simulation_function::Function;
                                 analysis_functions=Function[],
                                 storage_path="./results")
    
    storage_config = Dict(
        "path" => storage_path,
        "format" => "hdf5",
        "backup" => true,
        "compression" => true
    )
    
    workflow = ResearchWorkflow(name, description, parameter_space,
                               simulation_function, analysis_functions,
                               storage_config)
    
    # Create directory structure
    mkpath(storage_path)
    mkpath(joinpath(storage_path, "data"))
    mkpath(joinpath(storage_path, "figures"))
    mkpath(joinpath(storage_path, "reports"))
    
    # Save workflow definition
    workflow_file = joinpath(storage_path, "workflow_definition.json")
    save_workflow_definition(workflow, workflow_file)
    
    return workflow
end

"""
    save_workflow_definition(workflow::ResearchWorkflow, filename::String)

Save workflow definition for reproducibility.
"""
function save_workflow_definition(workflow::ResearchWorkflow, filename::String)
    definition = Dict(
        "name" => workflow.name,
        "description" => workflow.description,
        "timestamp" => string(now()),
        "julia_version" => string(VERSION),
        "parameter_space" => Dict(
            "parameters" => Dict(string(k) => string(v) for (k,v) in workflow.parameter_space.parameters),
            "fixed_params" => Dict(string(k) => v for (k,v) in workflow.parameter_space.fixed_params),
            "n_constraints" => length(workflow.parameter_space.constraints)
        ),
        "storage_config" => workflow.storage_config
    )
    
    open(filename, "w") do file
        JSON.print(file, definition, 2)
    end
end

# Demonstration examples
if abspath(PROGRAM_FILE) == @__FILE__
    println("🔬 Parameter Optimization - Research Examples")
    println("=" * 50)
    
    # Example 1: Simple parameter sweep for diffusion equation
    println("\n1. Diffusion Equation Parameter Study")
    
    # Define parameter space
    diffusion_params = Dict(
        :diffusion_coeff => 0.01:0.02:0.1,
        :grid_size => [32, 64, 128],
        :time_step => [0.001, 0.01, 0.1]
    )
    
    # Add stability constraint
    constraints = [
        p -> p.diffusion_coeff * p.time_step / (1.0/p.grid_size)^2 < 0.5  # CFL condition
    ]
    
    space1 = ParameterSpace(diffusion_params, constraints=constraints)
    
    # Simple simulation function (mock)
    function diffusion_simulation(params)
        # Mock simulation - returns convergence rate
        D, n, dt = params.diffusion_coeff, params.grid_size, params.time_step
        dx = 1.0 / n
        stability = D * dt / dx^2
        
        # Mock convergence analysis
        convergence_rate = stability < 0.25 ? 2.0 - 0.5*stability : 1.0 - stability
        error = dx^2 * (1 + 0.1*randn())  # Mock spatial error with noise
        
        return Dict(
            :convergence_rate => convergence_rate,
            :spatial_error => abs(error),
            :stability_number => stability,
            :computational_cost => n^2 * (1.0/dt)
        )
    end
    
    # Perform parameter sweep
    results1 = parameter_sweep(space1, diffusion_simulation, method=:grid, parallel=false)
    
    # Analyze results
    analysis1 = analyze_parameter_effects(results1, :spatial_error, plot_results=false)
    
    # Example 2: Random sampling for high-dimensional space
    println("\n2. High-Dimensional Random Sampling")
    
    # Larger parameter space
    high_dim_params = Dict(
        :param1 => 0.0:0.1:2.0,
        :param2 => 0.0:0.1:2.0,
        :param3 => 0.0:0.1:2.0,
        :param4 => 0.0:0.1:2.0,
        :param5 => 0.0:0.1:2.0
    )
    
    space2 = ParameterSpace(high_dim_params)
    
    # Mock optimization function
    function optimization_function(params)
        # Mock objective: minimize sum of squares with noise
        objective = sum(p^2 for p in values(params)) + 0.1*randn()
        return Dict(:objective => objective)
    end
    
    # Random sampling
    results2 = parameter_sweep(space2, optimization_function, 
                              method=:random, n_samples=200, parallel=false)
    
    analysis2 = analyze_parameter_effects(results2, :objective, plot_results=false, n_top=5)
    
    # Example 3: Latin hypercube sampling
    println("\n3. Latin Hypercube Sampling")
    
    results3 = parameter_sweep(space2, optimization_function,
                              method=:lhs, n_samples=100, parallel=false)
    
    analysis3 = analyze_parameter_effects(results3, :objective, plot_results=false, n_top=5)
    
    println("\n🎉 All parameter optimization examples completed!")
    println("Results summary:")
    println("  Grid search: $(length(results1.parameters)) combinations")
    println("  Random sampling: $(length(results2.parameters)) combinations") 
    println("  LHS sampling: $(length(results3.parameters)) combinations")
    println()
    println("Best objectives found:")
    println("  Random: $(round(minimum([r[:objective] for r in results2.results]), digits=6))")
    println("  LHS: $(round(minimum([r[:objective] for r in results3.results]), digits=6))")
    println("  (Theoretical minimum: 0.0)")
end